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Aerothermodynamic  Study  of  a  Generic  EFP  Configuration 


Kirk  J.  Vanden*  and  Douglas  V.  Nance1^ 
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The  present  study  addresses  the  aerodynamics  of  a  generic  explosively  formed  penetrator 
body  in  Mach  6  hypersonic  flight  at  sea  level.  As  a  first  study  of  this  type,  the  shape  of  the 
body  is  held  fixed,  and  the  flow  field  is  resolved  by  applying  state-of-the-art  large  eddy 
simulation  techniques  in  conjunction  with  a  hybrid  shock-turbulence  capturing  algorithm. 
Air  is  treated  as  a  mixture  of  nitrogen  and  oxygen,  and  the  governing  equations  are  closed 
by  a  modern  compressible  turbulence  closure  term.  Pressure  is  determined  by  using  the 
thermally  perfect  gas  equation  of  state  applied  to  each  species.  The  generic  body  also 
possesses  an  interior  cavity.  The  flow  field  is  captured  both  inside  and  outside  of  the  body. 
Since  aerodynamic  breakup  is  of  great  concern,  the  distribution  of  temperature  is 
determined  on  the  body  surfaces  as  well  as  temperature  gradients  based  upon  adiabatic  wall 
boundary  conditions.  Also,  the  structure  of  the  flow  field  is  examined  as  is  the  time  required 
for  stationarity.  These  factors  have  an  effect  on  flight  stability. 
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Nomenclature 

specific  internal  energy 
constant  pressure  specific  heat 
subgrid  kinetic  energy 
total  energy  per  unit  mass 
pressure 

heat  flux  vector  component 
species  k  perfect  gas  constant 
temperature 
time 

Cartesian  velocity  component 
Cartesian  coordinate 
delta  function 
density 

grid  scale  measure 
stress  tensor 


I.  Introduction 

AN  Explosively  Formed  Penetrator  (EFP)  is  composed  of  two  components,  an  explosive  lens  and  a  shaped  metal 
liner  as  is  shown  in  Figure  l.1  As  an  EFP  detonates,  the  detonation  wave  spreads  through  the  lens.  Controlled 
by  the  shape  of  the  lens,  the  wave  interacts  and  exerts  a  force  distribution  on  the  liner.  The  liner  accelerates  away 
from  the  blast,  but  the  center  region  travels  faster  than  the  remainder  of  the  liner.  As  a  result,  the  liner  takes  on  a 
crude  bell-like  shape.  The  distribution  of  explosive  force  causes  the  liner  to  undergo  large,  high  strain -rate 
deformation.  In  flight,  the  liner  drastically  changes  shape  from  a  disk  into  a  bell  shaped  bluff  body.  In  some  cases, 
an  open  cavity  may  exist  at  the  base  of  the  EFP.  Structurally,  the  liner  material  “flows”  in  the  plastic  regime,  and 
accordingly,  it  becomes  very  hot.  The  temperature  rise  due  to  the  effects  of  material  strain  is  complemented  by 
aerodynamic  heating  from  hypersonic  flow  in  the  thick  atmosphere  at  sea  level.  It  is  of  interest  to  investigate  the 
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Figure  1.  Cross-Section  for  a  notional 
Explosively  Formed  Penetrator. 


aerodynamics  of  a  generic  EFP  body  in  the  hypersonic  regime,  so  we  have  chosen  a  sea  level  velocity  corresponding 
to  Mach  6  since  this  Mach  number  is  comfortably  within  the  hypersonic  regime.  In  flight,  the  EFP  may  pitch  or  yaw 
in  an  unstable  way;  casting  a  complex  pattern  of  aerodynamic  loads  over  the  EFP  body.  These  forces  can  cause  the 
EFP  liner  to  break  apart  in  flight.  Investigating  the  phenomenon  of  aero  stability  and  aerodynamic  heating  are 
focuses  of  the  present  study. 

Simulating  the  events  comprising  EFP  formation  and  flight  require 
the  use  of  both  solid  and  fluid  mechanics  including  unsteady  hypersonic 
flow  physics  and  high  strain  rate  dynamics.  Our  current  capabilities  do 
„  not  permit  a  study  encompassing  the  structural  deformation  and  fluid 

£  ff  dynamics  together.  In  the  present  work,  the  structural  mechanics  of  the 

EFP  are  disregarded  while  we  concentrate  on  describing  the  flow  field 
around  a  canonical  EFP  body  configuration.  The  body  has  a  complex 
shape  characterized  by  a  blunt  nose  and  an  elongated  body  distending 
into  fin-shaped  lobes.  This  configuration  is  shown  in  Figure  2.  State-of- 
the-art  large  eddy  simulation  (FES)  techniques  are  used  to  simulate 
flow  around  the  body  for  a  fixed  Mach  number,  angle  of  attack  and 
angle  of  sideslip.  These  methods  allow  the  viscous,  turbulent  flow  field 
to  be  resolved  with  high  accuracy.  Moreover,  if  one  assumes  that  there 
is  no  heat  flow  through  the  thickness  of  the  liner,  the  temperature 
distribution  over  the  body  can  also  be  accurately  determined.  These 
flow  properties  can  be  determined  within  the  cavity  as  well  as  outside  of 
the  body  for  a  range  of  flight  conditions. 

During  the  period  of  time  encompassing  the  1960s  through  the 
1980s,  hypersonic  flow  was  intensively  studied  because  of  the  design  needs  of  reentry  vehicles  and  concept  designs 
such  as  the  National  Aerospace  Plane.2,3  In  deference  to  supersonic  flow,  there  is  no  strict  point  of  demarcation  in 
terms  of  speed  heralding  the  hypersonic  regime.  For  aerospace  vehicles  commonly  studied,  the  onset  of  hypersonic 
flow  physics  is  encountered  in  the  range  of  Mach  3  to  Mach  5,  or  even  higher.  The  physics  associated  with  this 

regime  is  very  different  that  that  associated  with 
slower  speed  flow.  Hypersonic  flow  is  characterized 
by  the  presence  of  thin  shock  layers  that  lie  between 
the  body  and  the  shock  wave.2  At  higher  Mach 
numbers,  this  layer  becomes  even  thinner,  and  it 
may  even  merge  with  the  boundary  layer.  The 
presence  of  an  entropy  layer  is  also  characteristic  of 
many  hypersonic  bluff  body  flow  fields.  Bluff  (or 
blunted)  bodies  cause  the  formation  of  curved  shock 
waves  that  stand  some  distance  away  from  the  body. 

The  strength  of  the  shock  wave  and  the  attendant 
entropy  jump  changes  along  with  shock  curvature.  It 
is  strongest  at  the  “nose”  and  weakens  as  the  shock 
wave  straightens  downstream.  This  effect  causes  a 
flow  with  strong  entropy  and  vorticity  gradients  to  wash  down  the  body  thickening  the  boundary  layer.2  This 
phenomenon  is  important  for  an  EFP  since  it  spends  much  of  its  flight  time  as  a  bluff  body  and  then  can  lengthen 
into  a  high  aspect  ratio  body  with  a  thin  shock  layer.  The  proximity  of  the  shock  wave  along  with  its  interaction  with 
the  boundary  layer  can  dictate  the  forces  and  moments  exerted  on  the  body.  These  considerations  have  a  significant 
effect  on  the  stability  of  the  EFP  in  flight. 

High  temperature  gas  dynamics  are  often  encountered  in  hypersonic  flow  fields.  In  order  to  satisfy  the  no -slip 
boundary  condition,  the  flow  must  reduce  from  hypersonic  speed  outside  of  the  boundary  layer  to  zero  at  the  body 
surface.  The  mechanism  decreasing  the  flow  speed  is  viscous  dissipation  or  friction.  The  kinetic  energy  of  the  flow 
is  transformed  into  the  internal  energy  of  gas  molecules  in  the  boundary  layer.  As  a  result,  boundary  layer 
temperature  increases  dramatically.2  In  slower  speed  flow  fields,  this  energy  is  absorbed  mostly  in  molecular 
translational  and  rotational  modes,  but  for  hypersonic  speeds,  energy  may  be  absorbed  in  vibrational  modes.  For 
sufficiently  high  speeds,  gas  molecules  may  either  ionize  or  dissociate  effectively  resulting  in  a  chemically  reacting 


Figure  2.  Geometric  configuration  for  a  generic  EFP. 
Forward  oblique  view  shown  on  the  left,  reverse  oblique 
view  shown  on  the  right. 
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flow  field.  The  gas  field  can  achieve  very  high  temperatures,  thousands  of  degrees  Kelvin.  This  thermal  energy  is 
transferred  in  part  to  the  EFP  body  softening  the  liner  material  and  creating  a  major  concern  for  weapon 
effectiveness. 

The  purpose  of  this  study  is  to  capture  some  of  these  effects  for  a  realistic,  yet  constant,  EFP  geometry.  Its 
velocity  is  held  at  constant  Mach  6  for  sea  level  flight  conditions.  It  is  in  this  sense  that  the  EFP  diverges  from 
common  hypersonic  flow  problems.  Both  re-entry  and  hypersonic  flight  occur  in  low  density  atmosphere.  At  sea 
level  density,  we  have  little  aerodynamic  data  for  this  speed  range.3  In  the  work  that  follows,  the  pitch  angle  is 
varied  as  a  parameter  over  the  range  of  zero  to  +7.5°.  For  each  pitch  angle,  the  time  required  to  achieve  a  stationary 
flow  field  is  determined.  The  existence  of  a  stationary  flow  indicates  that  flow  properties  vary  in  a  cyclical  manner. 
That  is  to  say,  regions  of  recirculation  may  exist  in  the  cavity  resulting  in  cyclical  pressure  fluctuations  on  the  body. 
As  a  result,  the  force  and  moment  distributions  alter  on  the  EFP  body.  The  flow  field  is  unsteady  but  varies  in  a 
periodic  manner.  The  settling  time  required  to  reach  this  state  is  also  of  theoretical  interest.  If  the  flow  remains  non¬ 
stationary  for  the  duration  of  the  flight,  the  body  may  continue  to  oscillate  in  pitch.  This  assertion  is  of  added 
significance  since  our  model  EFP  configuration  has  a  fixed  shape  for  the  duration  of  the  simulation.  The  temperature 
and  temperature  gradient  distributions  are  also  determined  on  the  interior  and  exterior  surfaces  of  the  body.  W e  also 
determine  whether  or  not  chemical  reactions  may  occur  in  the  flow  field  under  these  conditions.  Trends  in  the 
fluctuation  of  force  and  moment  are  calculated  at  the  different  pitch  angles. 


II.  Theory 

The  computer  program  Farge  Eddy  Simulation  with  FInear  Eddy  modeling  in  3  Dimension  (FESFIE3D)  has 
been  used  to  generate  the  EFP  flow  field  results  shown  later  in  the  paper.  FESFIE3D  is  developed  by  Suresh  Menon 
at  the  Georgia  Institute  of  Technology.4  This  computer  program  has  a  core  capability  for  dynamic  large  eddy 
simulation  (FES).  That  is  to  say,  flow  field  features  existing  at  the  scale  of  the  grid  or  larger  are  simulated  by 
numerical  solution  of  the  conservation  equations.  Flow  features  existing  at  the  subgrid  scales  (scales  smaller  than  the 
mesh  size)  must  be  modeled.5  The  subgrid  effects  are  represented  in  the  conservation  equations  via  closure  terms. 


A. Filtered  Governing  Equations 

The  conservation  laws  used  for  this  research  consist  of  the  compressible  Navier-Stokes  equations  cast  in  three 
dimensions.  For  real  turbulent  flow  fields,  we  cannot  solve  these  equations  via  direct  numerical  solution  due  to  the 
limitations  imposed  by  today’s  computer  resources.  Our  LES  approach  first  involves  spatially  filtering  the  governing 
equations  at  the  grid  scale  in  order  to  produce  a  set  of  turbulence  closure  terms  that  are  later  modeled.  The  filtered 
equations  are  written  as  follows.6 
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dt  dxt 
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where  k  =  1,. . Ns ,  Ns  being  the  number  of  gaseous  chemical  species  involved  in  the  calculation.  Equations  (1) 


through  (4)  are  mass,  momentum,  energy  and  species  conservation  equations,  respectively.  Symbols  p  ,  ui  ,  P  , 


Y k  and  qt  are  the  mixture  gas  density,  velocity  component  in  Cartesian  component  direction  x. ,  absolute 

pressure,  mass  fraction  for  the  kth  species  and  heat  flux  component  in  direction  xt ,  respectively.  Symbol  T 
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contains  the  components  of  the  Cartesian  shear  stress  tensor  while  E  represents  the  total  energy  per  unit  volume, 
i.e., 


Tij=V 


du du 


\dxi 


■  +  ■ 


dxt 

1 


,  du,  _ 
+  A — L&, 
dx, 


E  =  e  +  —  ui  ui 


(5) 


(6) 


The  summation  convention  is  applied  to  the  index  l  in  (5)  and  (6).  The  Fickian  diffusion  velocities  are  given  by 
Vt  k  .  The  overbar  A  notation  indicates  that  the  quantity  A  has  been  averaged  without  density  weighting  whereas 


the  notation  A  indicates  that  A  has  been  mass  averaged.  The  subgrid  stess  tensor  T-gs ,  subgrid  total  enthalpy  flux 

H*8S ,  subgrid  convective  species  flux  YtJs ,  subgrid  viscous  work  <J*8S  and  subgrid  species  diffusive  work  0*8k  are 
turbulence  closure  terms  that  require  modeling  to  some  degree.  Chemical  reactions  occurring  between  species  are 
governed  by  the  filtered  reaction  rate  term  (bk  .7  Pressure  in  the  gas  mixture  is  determined  through  the  use  of  a 
filtered  perfect  gas  equation  of  state,  i.e., 


P  =p(RT  +Tsgs) 


with  Rjj  equal  to  the  universal  gas  constant,  and 


(7) 


k= 1 


k= 1 


'  Ry  ' 

KMWkJ 


MWk  is  the  molecular  weight  of  species  k  .8  The  filtered  mixture  internal  energy  is  given  as 


s  s  1  ivs 


Sgs 
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k= 1  k= 1  o  k= 1 

In  equation  (9),  Cv  k  is  the  constant  volume  specific  heat  for  species  k  ,  and  Ek8S  is  the  subgrid  internal  energy  for 
species  k  .  The  heat  flux  term  may  be  written  as7 


ST  r  ~  ~ 

<li  =-K^  +  PLhkYkVi,k 

oxt  k=1 


where  K  is  the  local  coefficient  of  thermal  conductivity,  and8 


i 

K=Y,h%+Yt\CFl(J)dT 


(10) 


(11) 


The  specific  heat  at  constant  pressure  for  species  k  is  written  as 


Cp,k  ~  Cv,k  +  ^k 


(12) 


B.  Closure  Term  Modeling 

Taken  directly,  the  number  of  subgrid  terms,  denoted  by  superscript  SgS ,  result  in  a  system  of  governing 
equations  that  is  not  closed.  In  order  to  close  this  system  of  equations,  the  subgrid  terms  must  be  somehow  modeled. 
The  current  practice  is  to  neglect  Tsgs  ?  d-8k  and  Eskgs  since  these  terms  tend  to  be  small.6  The  crux  of  this 
modeling  effort  is  in  determining  the  subgrid  stress  tensor  based  upon  the  subgrid  kinetic  energy,  i.e., 
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rsgs  =  -2  vj  S  .  -  —Skk  8, ,  I  +  -ksgs8i , 

n  ij  3  kk  ij  i  3  ij 


where  the  evolution  equation  for  subgrid  kinetic  energy  is  given  by 
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The  subgrid  kinetic  energy  is  mathematically  defined  as 


sgs\3/2 


du.  (ksgs) 

Tir-^+pKe{-^4 


ksgs  =  —  (UjU,  -Ujilj )  (15) 

where  the  summation  convention  is  in  effect  for  the  index  l  ?  Equation  (12)  is  integrated  along  with  the 
conservation  equations,  and  since  it  is  model  equation,  it  relies  on  a  set  of  parameters  such  as  Vt ,  the  eddy  viscosity, 
given  by 

v,  =K„VT^A  (16) 

where  Kv  is  a  coefficient  set  by  the  Locally  Dynamic  subgrid  Kinetic  energy  Model  (LD KM). 8,10  The  dissipation 
term  for  ksgs  is  given  by 

^  pK(ksgs)312 

D  =H  gV  (17) 

K  A 

with  as  the  LDKM  dissipation  coefficient  (analogous  to  );  A  is  the  local  grid  scale  measurement.  Other 
parameters  introduced  in  (12)  are  jU  the  dynamic  viscosity,  Pr?  the  turbulent  Prandtl  number,  Msgs  the  turbulent 
Mach  number  based  upon  ksgs ,  and  GC  d  is  a  pressure-dilatational  scaling  coefficient.  Note  also  that  (13)  and  (14) 

are  tightly  coupled  through  the  presence  of  both  ksgs  and  T-gs . 

The  theory  that  underlies  LDKM  is  quite  complex,  and  the  details  of  implementing  the  above  equations  (along 
with  others  not  shown)  are  beyond  the  scope  of  this  paper.  However,  it  is  worthwhile  to  discuss  the  physical 
meaning  behind  this  model.  In  the  first  place,  one  may  note  the  number  of  coefficients  existing  in  (11),  (12),  (14) 
and  (15).  The  proliferation  of  coefficients  indicates  that  LDKM  is  indeed  a  model,  but  these  coefficients  are  not 
tunable.  Instead,  they  are  in  most  cases  based  upon  knowledge  of  the  dynamic  behavior  of  turbulent  compressible 
flow  fields.  That  is  to  say,  these  coefficients  are  set  either  autonomously  by  LDKM  operating  in  dynamic  mode,  or 
they  have  fixed  or  easily  calculated  values.10  Also,  LDKM  diverges  from  most  contemporary  dynamic  LES  models 
through  the  presence  of  the  term  marked  by  an  (*)  asterisk.  This  term  represents  dilatational  effects  within  the  flow 
field  due  to  pressure.  It  is  through  this  term  that  (14)  truly  becomes  representative  of  highly  compressible  flow 
fields,  so  this  equation  stands  as  a  hallmark  for  state-of-the-art  research  in  this  field.  When  operating  LDKM  in 
dynamic  mode,  the  subgrid  scale  properties  are  made  to  behave  in  an  inertially  correct  manner  with  respect  to  the 
resolved  flow  field. 
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III.  Numerical  Methods 

In  the  preceding  section,  the  Navier-Stokes  equations  are  presented  along  with  the  equations  for  Menon’s 
LDKM.  These  equations  are  well  suited  for  describing  turbulent  flow  fields  in  any  flight  regime  (subsonic, 
supersonic  and  hypersonic).  In  this  section,  we  discuss  the  numerical  algorithms  associated  with  solving  these 
equations.  The  EFP  problem  requires  hypersonic  flow  solutions  characterized  by  the  presence  of  strong  shock  waves 
and  high  temperature  gas  dynamics.  It  is  difficult  for  computer  codes  to  at  the  same  time  capture  both  shock  waves 
and  subtle  turbulent  flow  field  fluctuations  because  each  of  these  phenomena  requires  the  use  of  different  numerical 
algorithms.  Generally,  numerical  dissipation  must  be  used  to  dampen  oscillations  around  shock  waves,  yet  this  same 
dissipation  washed  away  the  weaker  vortical  motions  associated  with  turbulence.  Moreover,  high  order  centered 
difference  stencils  are  usually  required  to  accurately  capture  turbulence,  but  these  stencils  cause  severe  spurious 
oscillations  around  shock  waves.  It  is  also  important  to  realize  that  the  shock-capturing  algorithms  must  be  able  to 
contend  with  complex  equations  of  state.  You  may  recall  that  we  must  apply  the  thermally  perfect  gas  equation  of 
state  in  the  hypersonic  flight  regime.  After  performing  several  years  of  research,  we  have  arrived  at  a  numerical 
scheme  that  is  eclectically  composed  of  two  parts.  For  shock-capturing,  we  apply  what  is  referred  to  as  the  Harten- 
Fax-van  Feer-Contact/Einfeldt  (HFFC/E)  approximate  Riemann  solver.11,12  In  smooth  regions  of  the  flow  field,  we 
apply  a  second/fourth  order  MacCormack  solver  used  in  the  earlier  incarnations  of  FESFIE3D.  A  switching 
algorithm  allows  the  computer  program  to  autonomously  select  the  appropriate  space  scheme  based  upon  local 
properties  throughout  the  flow  field.  Below,  we  present  brief  descriptions  of  the  spatial  integration  schemes. 

A.  The  HLLC/E  Approximate  Riemann  Solver 

HEEC/E  is  composed  of  HLLE  (HLL-Einfeldt),  a  scheme  that  does  not  preserve  the  contact  wave,  and  HLLC 
(HLL-Contact),  a  scheme  that  does  capture  the  contact  wave.  In  this  context,  the  contact  wave  is,  in  fact,  a  contact 
discontinuity,  a  common  flow  feature  exhibited  by  the  shock  tube  problem  as  well  as  by  strong  blast  waves.  These 
Godunov-like  schemes  are  constructed  on  a  one-dimensional  coordinate  oriented  normal  to  a  finite  volume  interface. 
It  is  in  this  manner  that  these  schemes  are  applicable  to  three-dimensional  flow  fields.  The  numerical  solution  is 
calculated  for  a  vector  of  conserved  variables 

U  =  {p,u,v,w,pE,pksgs,pYv...,pYN)  (18) 

the  associated  flux  vector 


F  =  (pq,puq  +  Pnx , pvq  +  Pn  , pwq  +  Pnz , ( pE  +  P)q, pksgsq, pY^q,.. pYN  q) 


(19) 


where  q  =  unx  +  vny  +  wnz ,  and  h=nxi  +  ny  j  +  n7k  is  the  unit  surface  normal  vector  at  the  interface.  Based 
upon  these  definitions,  the  HLLE  numerical  flux  vector  can  be  written  as13 


r'HLLE 

A+1/2 


Fl,  0  <SL 

F\  SL<0<SR 
Fr,  Sr<  0 


(20) 


where 

F*  = 

The  Einfeldt  wave  speed  estimates  are 


SRFL -SLFR  +SLSR(UR - UL ) 


SR -SL 


(21) 


SL  =min(^L  —  cL,q  —  c)\  SR  =max(qR  +  cR  ,q +  c)  (22) 

In  (22),  c  is  the  speed  of  sound,  and  the  notation  ( w )  indicates  the  Roe  average.14  This  numerical  flux  works  quite 
well  for  hypersonic  flow  and  is  not  susceptible  to  contact  instabilities  that  can  occur  when  strong  bow  shocks  are 
aligned  with  the  grid.  Unfortunately,  because  of  dissipative  effects  near  the  contact  wave,  it  cannot  be  used  for 
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extended  regions  of  turbulent  flow.13  In  turbulent  regions  of  the  flow  field,  the  HLLC  numerical  flux  shows 
improved  performance.  The  numerical  flux  for  the  HLLC  scheme  has  a  different  form,  i.e., 

Fl  ,  0<SL 

FL*  =FL  +SL(UL* -UL),  SL  <0<S* 

\Fr*  =Fr  +  Sr(Ur* -Ur),  sl*<o<sr 
Fr  ,  SR  <0 


FHLLC  _ 
ri+M2  -  s 


(23) 


In  (23),  intermediate  states  U  and  U  have  been  introduced  in  order  to  model  the  presence  of  the  contact 
discontinuity.  These  states  may  be  written  as13 


UL*  =aLUL  +  (0,pLa>Lnx,pLa>Ln  ,pLa>Lnv,i//L,0,...) 


with 


«L=fiL+‘l-  aL=-PL(qL-SL) 


(24) 


(25) 


and  the  intermediate  wave  speed 

s*  =  ■ 


PL  .  L  L/cL  L\  K  K/cK  A\ 

+£  q  (5  -q  )-p  q_(S_-q_) 

^L/riL  L\  R  /  CR  R\ 

p  ( S  -q  )-p  ( S  -q  ) 


(26) 


The  state  UR  is  formed  by  replacing  R  for  L  in  (24)  through  (26).  In  LESLIE3D,  both  of  these  numerical  flux 
formulations  are  used  to  form  the  HLLC/E  scheme.  We  state  the  result  only,  i.e., 


^  HLLC  IE 

A+1/2 


j  Clv2  >  ( dpj  <  0  and  du  j  <  0)  or  (dP  k  <  0  and  du  k  <  0) 

{L,+v2  ’  otherwise 


(27) 


This  method  requires  shock  detection  in  directions  transverse  to  the  direction  of  the  interface  normal,  so 

y 

,.^-u)  3 


j  _  I  Pj,j+U  ^jpjc  I 

PJ  min  (P 


i,j+U  ’ 


(28) 


du,j  ~  Ui,j+'\,k  Ui,j--\,k 


(29) 


and  so  forth. 

A  close  inspection  of  the  numerical  flux  formulas  above  reveals  that  the  computation  of  the  numerical  flux  first 
requires  the  left  ( L )  and  right  ( R  )  states  be  identified  at  interface  location  ( /  + 1  /2).  These  “upwind”  states  are 
reconstructed  from  the  date  stored  in  the  finite  volume  cells.  This  process  is  accomplished  with  the  use  of 
Monotone  Upwind  Scheme  for  Conservation  Laws  (MUSCL)  interpolation.14  High  order  interpolation  also  requires 
the  use  of  a  nonlinear  limiter  with  “flattening”15  in  order  to  maintain  data  monotonicity.13  The  details  of  these 
procedures  are  omitted  from  this  work;  interested  readers  are  initially  referred  to  References  13  through  15. 


B.  The  MacCormack  Solver 

For  the  “hybrid”  solver,  that  is  suitable  for  capturing  shock  waves  and  turbulence  anywhere  in  the  flow  field,  it  is 
necessary  to  switch  between  the  upwind  scheme  and  a  centralized  scheme  based  upon  MacCormack’ s  method  when 
calculating  the  numerical  flux  at  interfaces.  The  finite  difference  version  of  this  scheme  uses  forward  and  backward 
differences  alternatively  to  remove  bias  from  numerical  error.  For  the  finite  volume  scheme,  we  choose  forward  and 
backward  upwind  variables  alternatively  for  the  flux  computation,  i.e., 

(30) 


pNum  =  J7(TJ±  ) 

1  i+M2  1  \ui+M2S 


where 
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U  mi  2  =  ;  Ui+ 1/2  =  t/£-  (31) 

It  can  be  shown  that  this  scheme  is  stable  and  second  order  accurate  in  space.13  LESLIE3D  also  offers  a 
MacCormack  scheme  accurate  to  the  fourth  order.  The  hybrid  method  is  required  to  switch  back  and  forth  between 
the  upwind  and  MacCormack  solvers.  The  “switch”  is  based  upon  the  “smoothness”  of  the  flow  field.8  Principally, 
LESLIE3D  autonomously  switches  on  the  upwind  scheme  only  in  the  vicinity  of  flow  field  discontinuities  and 
employs  the  MacCormack  solver  in  smooth  regions  of  the  flow  field.  This  hybrid  methodology  performs  quite  well 
and  has  been  extensively  validated  for  the  Richtmyer-Meshkov  instability,  the  chosen  test  problem  for  this  research.9 
It  has  shown  a  great  deal  of  efficacy  for  simulating  shock-turbulence  interaction  problems.7 


C.  Time  Integration  Scheme  and  Computer  Code  Structure 

LESLIE3D  utilizes  explicit  time  integration  to  capture  the  physics  associated  with  unsteady  flow  fields. 
Although,  the  explicit  time  step  is  limited  by  the  Courant-Friedrich-Lewy  (CFL)  criteria,  it  presents  a  high  level  of 
efficiency  and  computational  simplicity  on  parallel  machines.  A  version  of  the  modified  Euler  method  is  applied  for 
the  MacCormack,  upwind  and  hybrid  schemes.  This  time  integration  method  (in  two  steps)  is  succinctly  summarized 
for  one  space  dimension  as  follows.16 


U*=  t/"_1+  — 
Ax 


Num.n-A 

+1/2 


jjNum,n- 1  ' 


i-M  2 


(32) 


Un 


u  +u 


n- 1 


At 

Ax 


(C+1/2 


T?Num*  \ 
^  i-M  2  ) 


(33) 


Equations  (32)  and  (33)  require  minor  adaptations  for  use  with  problems  cast  in  three  dimensions.  Basically,  the 
length  term  Ax  roughly  corresponds  to  a  ratio  of  the  cell  volume  to  an  average  surface  area  for  the  interface  normal 
to  the  coordinate  indexed  by  i .  The  advective  terms  in  (1)  through  (4)  are  easily  discretized  for  structured, 
hexahedral  finite  volumes  while  the  viscous  terms  require  a  basic  generalized  coordinate  transformation.  The 
mathematical  procedures  are  relatively  easy  to  accomplish  and  are  quite  robust. 

It  is  prudent  to  discuss  the  larger  aspects  of  LESLIE3D  as  well.  LESLIE3D  is  a  finite  volume  solver  written  for 
arbitrary  structured  grids  and  is  fully  parallelized.  The  current  version  of  the  code  is  implemented  in  multi -block 
form  for  analyzing  highly  complex  flow  configurations.  LESLIE3D  is  also  equipped  with  finite  rate  chemistry 
algorithms  with  a  choice  of  turbulent  chemical  closure  schemes.  For  problems  involving  evaporating  reactive 
droplets  or  solid  particles,  LESLIE3D  offers  a  dispersed  phase  that  is  fully  coupled  with  the  gas  phase.  A  dense 
dispersed  phase  algorithm  has  recently  been  implemented  within  the  computer  program  lending  to  its  capability  for 
solving  multiphase  physics  problems. 


IV.  TEST  PROBLEMS  AND  RESULTS 


A.  Problem  Set-up  and  Basic  Aerostability  Analyses 

A  hypersonic  flow  field  consisting  of  nitrogen  and  oxygen  surrounding  the  generic  EFP  configuration  (shown  in 
Figure  2)  has  been  analyzed.  In  this  example,  the  EFP  travels  at  Mach  6  under  normal  atmospheric  conditions  at  sea 
level.  Four  pitch  angles  have  been  analyzed  for  this  study:  0°,  2.5°,  5°  and  7.5°.  Other  flight  attitudes  are  certainly 
possible  (including  yaw  angles),  but  the  constraints  of  computer  availability  have  only  permitted  the  time  required  to 
consider  these  four  conditions.  The  lobed  design  of  the  EFP  complicates  the  flow  field  geometry.  In  addition  to  the 
exterior  region,  the  flow  field  is  also  computed  inside  of  the  EFP  folds;  see  Figure  2.  The  EFP  body  is  non- 
deforming,  and  the  pitch  angle  far  field  conditions  are  implemented  along  the  outer  boundary  of  the  grid.  Since  the 
freestream  is  supersonic,  the  enforcement  of  these  conditions  presents  no  difficulty  for  the  range  of  pitch  angles 
under  consideration.  The  grid  is  broken  into  243  subdomains  and  consists  of  just  over  4.6  million  grid  points.  For 
efficient  parallel  computing,  the  subdomains  are  balanced  to  operate  on  144  processors.  For  each  pitch  case,  the 
flow  field  is  “impulsively”  started  by  initializing  the  flow  in  each  cell  to  a  velocity  of  1 00  m/s  in  the  direction  of  the 
freestream.  As  LESLIE3D  operates,  the  flow  field  then  stabilizes  to  a  stationary  solution  around  the  EFB  body 
permitting  the  collection  of  physical  information  including  the  temperature  distribution  on  the  body. 
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Figure  3.  Components  of  total  force  in  Newtons  exerted  on  the  EFP  at  (left)  0°  pitch  and  (right)  2.5°  pitch. 
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Figure  4.  Components  of  total  force  in  Newtons  exerted  on  the  EFP  at  (left)  5°  pitch  and  (right)  7.5°  pitch. 


The  settling  time  for  the  flow  field  is  of  great  interest  since  time  is  required  for  the  flow  field  to  stabilize  around 
a  real  EFP  that  has  just  been  fired.  Figure  3  contains  plots  of  the  total  force  components  exerted  on  the  EFP  in  the  X 
for  pitch  angles  0°  and  2.5°,  respectively.  Figure  4  contains  corresponding  plots  for  the  5°  and  7.5°  pitch  cases.  As 
one  can  see,  for  the  5°  and  7.5°  cases  (Figure  4),  the  flow  field  does  not  begin  to  exhibit  stationarity  until  over  1  ms 
of  flight  time  has  elapsed.  The  same  trend  is  exhibited  in  Figure  3,  but  as  this  manuscript  is  being  drafted,  the 
numerical  solutions  are  still  in  progress  for  the  0°  and  2.5°  pitch  angles.  At  Mach  6,  the  elapsed  flight  time  of  one 
millisecond  corresponds  to  about  28  body  lengths  or  2.25  meters  of  travel.  For  earlier  times,  the  flow  field  around 
the  body  is  expected  to  be  highly  unsteady,  but  non-inertial  calculations  are  needed  to  prove  this  assertion.  In  Figure 
5,  total  force  exerted  in  the  X  direction  is  shown  for  each  of  the  pitch  cases.  For  the  7.5°  case,  the  force  exceeds  that 
exhibited  at  5°  pitch  due  to  the  increased  angle  of  attack,  i.e.,  more  of  the  body’s  surface  area  is  exposed  to  increased 
pressure  drag.  It  is  expected  that  this  trend  will  also  be  reflected  in  the  0°  and  2.5°  pitch  cases.  However,  due  to  the 
restrictions  of  time,  the  computer  solution  for  the  2.5°  case  has  not  yet  advanced  far  enough  to  prove  this  assertion. 
Still,  from  the  trend  observed  thus  far,  it  seems  unlikely  that  a  conventional  EFP  will  achieve  stationary  flow 
conditions  within  its  normal  time  of  flight.  If  this  assertion  bears  true,  it  may  have  interesting  implications  for 
aerostability.  The  inherent  unsteadiness  of  the  flow  field  may  complicate  the  process  of  stabilizing  the  EFP’s  flight. 
We  can  also  differentiate  between  the  effects  of  pressure  and  viscous  forces.  The  generic  EFP  geometry  is  still  a 
streamlined  body,  so  it  is  expected  that  the  force  due  to  pressure  exceeds  that  of  skin  friction  drag.  Generally,  we  do 
not  speak  of  “lift”  for  this  type  of  body,  but  both  drag  and  significant  side  forces  do  arise.  (See  Figures  3  and  4). 
This  assertion  is  illustrated  by  Figure  6,  a  plot  comparing  total  force  with  skin  friction  drag  (or  viscous  force)  for 
pitch  angles  0°,  2.5°,  5°  and  7.5°.  From  this  graph,  we  can  conclude  that  the  forces  exerted  by  pressure  (including  the 
shock  wave  drag)  overwhelm  viscous  forces  exerted  on  the  EFP.  The  total  force  magnitude  increases  along  with  the 
pitch  angle.  Interestingly,  the  viscous  force  retains  essentially  the  same  magnitude  in  time  as  the  pitch  angles 
changes.  It  follows  that  the  fluctuating  pressure  force  dominates  the  motion  of  the  EFP,  at  least  for  the  lower  pitch 
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angles.  Given  that  pressure  forces  are  usually  higher  than  skin  friction  drag,  it  is  expected  that  the  same  result  will 
still  apply  at  higher  pitch  angles. 


Time  (ms) 

Figure  5.  X-Component  of  aerodynamic  force  in  Newtons  exerted  on  the  EFP  at  all  pitch  angles. 
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Figure  6.  Total  force  magnitude  in  Newtons  compared  with  viscous  force  magnitude  for  all  pitch  angles. 

B.  Turbulent  Flow  Field  Physics 

Typical  hypersonic  flow  characteristics  are  exhibited  near  the  blunt  nosed  EFP.  The  proximity  of  the  detached 
bow  shock  wave  motivates  an  examination  of  the  flow  field  around  the  body.  Our  hybrid  LES  simulation  easily 
affords  doing  so,  but  it  is  important  to  bear  in  mind  the  limitations  of  3D  field  post-processing.  Recall  that  the 
computational  grid  consists  of  243  subdomains;  each  subdomain  consists  of  one  structured  grid,  but  there  is  no 
global  index  system,  (i.e.,  an  (/ ,  j  ,k)  ordered  triple)  available  to  address  the  grid  on  the  whole.  To  accurately 
extract  data  from  the  flow  field,  we  normally  plot  along  constant  index  planes  to  obtain  a  smooth  flow  field  property 
distribution.  For  complex  domain  decompositions,  it  is  not  possible  to  plot  data  smoothly  along  such  planes.  As  a 
result,  the  post-processor  is  required  to  interpolate  between  constant  index  planes  over  more  than  one  subdomain. 
This  deficiency  can  cause  the  plotting  program  to  create  uneven  contours.  This  difficulty  is  illustrated  in  Figures  7 
and  8,  plots  of  flow  field  properties  generated  at  the  nose  of  the  EFP.  Figure  7  contain  plots  of  gas  density  and 
pressure.  Although  the  contours  in  these  figures  suffer  from  gaps  due  to  interpolation,  we  can  clearly  see  the 
formation  of  a  bow  shock  that  stands  away  from  the  nose.  The  shock  wave  is  identified  as  a  yellow-red  high  density 
region  separated  from  the  nose  by  a  blue  low  density  region.  Moreover,  Figure  7  also  shows  the  high  pressure  post¬ 
shock  region  very  clearly.  Figure  8  illustrates  the  region  of  high  vorticity  behind  the  curved  shock,  and  we  can  also 
see  subgrid  kinetic  energy  spike  in  this  region.  The  results  conveyed  by  Figure  8  are  very  important;  they  clearly 
convey  the  presence  of  strong  turbulent  eddies  near  the  EFP  body.  Both  vorticity  and  subgrid  kinetic  energy  take  on 
a  wide  range  of  values  at  the  nose.  The  resulting  strong  gradients  in  these  quantities  advect  along  the  surface  of  the 
EFP.  Hence,  the  EFP  boundary  layer  is  turbulent.  If  turbulence  were  not  generated  at  the  nose,  we  would  see 
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Figure  7.  Contours  of  flow  field  (left)  density  in  kg/m3  and  (right)  pressure  in  Pascals  generated  along  the 
central  z  plane  at  the  nose  of  the  EFP.  Plots  correspond  to  1.35  ms  of  elapsed  time  for  the  5°  pitch  case. 


Figure  8.  Contours  of  flow  field  (left)  vorticity  in  s'1  and  (right)  subgrid  kinetic  energy  in  m2/s2  cast  along  the 
central  z  plane  at  the  nose  of  the  EFP.  Plots  correspond  to  1.35  ms  of  elapsed  time  for  the  5°  pitch  case. 

turbulence  die  out  as  the  solution  progresses.  This  behavior  would  be  assured  by  the  dissipation  term  in  Equation 
(14).  Figure  8  is  especially  significant  when  viewed  along  with  Figures  9  and  10.  Figure  9  contains  plots  of  subgrid 
kinetic  energy  taken  at  the  first  cell  adjacent  to  the  external  body  surface.  The  first  plot  is  generated  along  a 
longitudinal  curve  starting  at  the  nose  of  the  EFP  and  terminating  at  its  trailing  end.  The  curve  is  confined  to  the 
XZ  -plane  on  the  “top”  surface  of  the  EFP.  The  second  plot  is  generated  along  an  azimuthal  ring  circling  the  EFP 
body  parallel  to  the  yz  -plane  at  the  location  X  =3.659  mm.  The  existence  of  turbulence  is  clearly  indicated  in 
Figure  9.  The  longitude  plot  shows  strong  fluctuations  in  subgrid  kinetic  energy  along  the  length  of  the  EFP  body 
with  the  strongest  turbulent  eddies  existing  behind  the  nose  and  near  the  “shoulder  region”  for  the  EFP.  The 
magnitude  of  the  turbulence  varies  by  as  much  as  a  factor  of  ten.  Between  the  two  solutions  shown  in  time,  one  can 
see  a  disturbance  in  subgrid  kinetic  energy  propagating  from  the  tail  toward  the  nose  in  the  subsonic  boundary  layer. 
The  azimuth  plot  (left  in  Figure  9)  shows  significant  changes  in  this  energy  around  the  body  occurring  over  time. 
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Longitude  Station  Azimuthal  Station 

Figure  9.  Subgrid  kinetic  energy  (J/kg)  calculated  in  the  first  cell  adjacent  to  the  body:  (left)  along  a 
longitudinal  curve  from  nose  to  trailing  end,  and  (right)  along  an  azimuthal  ring  a  short  distance  behind  the 
nose  for  the  5°  pitch  case. 


Longitude  Station  Azimuthal  Station 

Figure  10.  Vorticity  magnitude  (s'1)  calculated  in  the  first  cell  adjacent  to  the  body:  (left)  along  a  longitudinal 
curve  from  nose  to  trailing  end,  and  (right)  along  an  azimuthal  ring  behind  the  nose  for  the  5°  pitch  case. 


Interestingly,  the  strongest  fluctuations  occur  on  the  “sides”  of  the  EFP  subjected  to  cross-flow  conditions.  Based 
upon  the  behavior  of  subgrid  kinetic  energy  in  these  regions,  it  is  worthwhile  to  examine  the  magnitude  of  vorticity 
shown  in  Figure  10.  The  range  of  vorticity  magnitude  is  so  wide,  we  have  employed  a  log  scale  for  plotting.  Note 
that  the  longitude  plot  also  reveals  the  existence  of  a  disturbance  propagating  upstream  in  the  boundary  layer. 
Vorticity  also  spikes  in  essentially  the  same  locations  predicted  by  Figure  9.  The  azimuth  plot  indicates  that  in  time 
vorticity  fluctuates  strongly  around  the  body  especially  in  areas  with  strong  fluctuations  in  subgrid  turbulent  kinetic 
energy.  By  revisiting  Figure  8  (left),  we  observe  that  bubbles  of  high  vorticity  develop  in  the  node  and  shoulder 
regions  and  then  advect  along  the  body  surface.  This  information  confirms  that  a  turbulent  flow  field  driven  by 
localized  sources  exists  around  the  EFP  during  its  flight.  Also,  vorticity  seems  to  be  coupled  with  gradients  in 
subgrid  kinetic  energy.  For  these  reasons,  it  is  doubtful  that  the  EFP  flow  field  can  attain  stationary  flow  conditions 
within  its  short  time  of  flight. 


C.  The  Distribution  of  Temperature 

Hypersonic  bluff  body  flow  fields  are  characterized  by  the  presence  of  high  temperature  gas  dynamics.  Slices  of 
the  temperature  flow  field  are  presented  Figure  1 1  for  the  5°  pitch  case.  The  presence  of  a  high  temperature  bubble  is 
evident  at  the  EFP  nose.  The  temperature  in  this  neighborhood  rises  above  2000°K;  as  a  result,  the  gases  become  hot 
enough  to  support  limited  chemical  ionization.3  This  assertion  is  also  supported  by  Figure  8  since  the  subgrid  kinetic 
energy,  a  driver  for  chemical  reactivity,  becomes  very  high  in  this  region.  The  strong 
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Figure  11.  Flow  field  temperature  in  Kelvin  along  the  central  vertical  mid-plane  at  (left)  the  nose  and  (right) 
for  the  nose  and  shoulder  region.  In  the  plot  on  the  right,  the  large  blue  region  is  for  the  interior  of  the  EFP. 
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Figure  12.  EFP  Surface  Area  Fraction  versus  surface 
temperature  (K)  for  all  pitch  angles. 


correlation  between  temperature  and  subgrid  kinetic 
energy  lends  to  the  likelihood  of  chemical  reactions. 
When  time  permits,  this  problem  may  be  solved  using 
turbulent  chemistry  algorithms  since  the  attendant 
changes  to  the  flow  field  can  affect  EFP  motion.  It  is 
interesting  to  analytically  examine  the  distribution  of 
temperature  over  the  body  surface  at  different  pitch 
angles.  Figure  12  addresses  this  subject  by  dividing  the 
temperature  range  from  200 °K  to  3000  °K  into  bins.  We 
then  sum  the  incremental  areas  along  the  EFP  surface 
that  are  “wetted”  by  temperature  within  the  respective 
bins.  This  Figure  illustrates  the  fraction  of  EFP  surface 
area  subjected  to  temperatures  within  the  overall  range. 
Although  the  solution  corresponding  to  the  2.5°  pitch 
angle  has  not  fully  matured,  Figure  12  reveals  useful 
information.  We  note  that  there  is  a  strong 
concentration  of  surface  area  corresponding  to  the 
temperature  band  between  2000°K  and  2100°K.  In  this 


range,  ionization  reactions  become  likely  for  Oxygen.  As  a  result,  this  calculation  is  to  be  performed  with  finite  rate 
turbulent  chemistry  algorithms  activated.  It  is  possible  that  the  surface  of  the  EFP  may  also  be  wetted  by  an  ion 
field.  If  so,  the  flow  properties  around  the  body  would  be  altered.  In  any  case,  the  higher  gas  temperatures  indicated 
near  the  body  motivate  a  detailed  examination  of  temperature  along  the  body  surface.  An  EFP’s  in-flight 
survivability  is  largely  a  function  of  its  surface  temperature  since  the  structural  integrity  of  the  EFP  depends  on 
thermal  softening  of  the  metal  liner.  If  regions  of  high  temperature  occur  on  the  body,  the  EFP  may,  under 
aerodynamic  loading,  break  apart  in  flight.  For  the  5  °  pitch  case,  the  external  and  internal  surface  temperature  plots 
are  shown  in  Figure  13  at  1.21  ms  solution  time.  Similar  behavior  is  observed  at  other  pitch  angles.  As  one  may 
expect,  there  are  large  differences  in  temperature  between  the  exterior  and  interior  surfaces.  It  is  interesting  to  see 
that  the  high  and  low  temperature  spatial  distributions  nearly  reverse  between  the  exterior  and  interior.  This 
phenomenon  is  likely  to  be  due  to  recirculating  flow  at  the  base  of  the  EFP.  Of  course,  the  interior  exists  at  cooler 
temperatures  that  the  exterior  surface  since  our  model  does  not  capture  conductive  heat  transfer  through  the  metal 
liner.  Still,  we  can  easily  calculate  temperature  gradients  on  these  surfaces.  The  plots  are  shown  in  Figure  14.  The 
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Figure  13.  Temperature  (K)  distributions  on  the  (left)  exterior  and  (right)  interior  surfaces  of  the  EFP  at  5° 
pitch  angle,  solution  time  1.21  ms. 


Figure  14.  Magnitude  of  temperature  gradient  (K/m)  distributions  on  the  (left)  exterior  and  (right)  interior 
surfaces  of  the  EFP  at  5°  pitch  angle,  solution  time  1.21  ms. 

magnitude  of  the  temperature  gradient  on  the  exterior  surface  is,  to  say  the  very  least,  intense,  particularly  near  the 
oblique  surfaces  of  the  nose,  on  the  shoulders  and  on  the  longitudinal  ribs  of  the  EFP.  The  differential  change  in  the 
gradient  is  less  pronounced  on  the  interior  surface,  but  the  greater  part  of  the  EFP  body  surface  is  subjected  to  very 
high  thermal  gradients  with  associated  thermal  stress.  It  is  also  important  to  realize  that  any  irregularity  (or  crease) 
in  or  on  the  body  results  in  a  region  of  high  temperature  due  to  frictional  heating.  Note  that  in  Figure  13  there  is 
ring-shaped  defect  just  behind  the  nose  that  is  subject  to  temperatures  exceeding  2200°K.  The  fact  is  important  since 
for  real  EFPs,  body  surface  irregularities  are  commonplace.  The  strong  localization  of  temperature  indicated  leads 
one  to  wonder  at  how  other  properties  might  be  distributed  over  the  EFP  surface.  Consider  vorticity  magnitude  as  an 
example.  This  property  is  presented  in  Figure  15  for  the  exterior  surface.  As  is  indicated  by  this  plot,  vorticity 
changes  by  over  a  factor  of  10  on  the  body.  Hence,  the  boundary  layer  is  subject  to  violent  mixing,  and  the  spotty 


14 


Approved  for  Public  Release;  Distribution  Unlimited  (96ABW-2009-0376;  Sept  8, 2009) 


Figure  15.  Distribution  of  vorticity  magnitude  (s-1)  on  the  exterior  EFP  surfaces  at  5°  pitch  angle,  solution 
time  1.21  ms. 

turbulent  fluctuations  in  vorticity  are  clearly  evident.  Strong  vorticity  fluctuations  can  affect  the  pressure  distribution 
at  the  body  surface  and  impact  aerodynamic  loading.  Hence,  predicted  EFP  motion  can  be  altered  by  the  accuracy  of 
the  aerodynamic  model.  An  important  conclusion  is  that,  due  to  this  mixing,  accurate  turbulence  models  must  be 
used  in  order  to  accurately  predict  the  temperature  distribution  for  the  EFP  in  time. 


IV.  CONCLUSIONS 

In  the  course  of  this  research,  we  have  performed  a  time-accurate  simulation  of  the  time-dependent  hypersonic, 
impulsively  started  flow  field  evolving  around  a  generic  EFP  shape.  The  example  simulation,  performed  with  the 
use  of  the  LESLIE3D  multiphase  physics  computer  code,  has  been  conducted  at  Mach  6  in  a  normal  atmospheric 
mixture  of  Oxygen  and  Nitrogen  under  sea  level  flight  conditions  for  a  series  of  pitch  angles.  LESLIE3D  is  a 
massively  parallel  computer  code  capable  of  the  high  fidelity  resolution  of  shock-turbulence  including  chemical 
reactions.  Our  goal  has  been  to  capture  the  physics  associated  with  this  flow  field  as  accurately  as  possible  by  using 
state-of-the-art  LES  techniques.  The  temperature  field  is  resolved  with  the  use  of  the  Thermally  Perfect  Gas 
equation  of  state.  A  principal  conclusion  of  this  research  is  that  the  EFP  flow  field  is  characterized  by  the  presence 
of  a  highly  turbulent  boundary  layer.  Effects  of  this  turbulence  must  be  accounted  for  in  the  calculation  of  loads  due 
to  the  interaction  of  subgrid  kinetic  energy  with  regions  of  high  vorticity.  Fluctuating  fields  of  subgrid  kinetic 
energy  advect  through  the  boundary  layer;  upstream  feedback  is  also  noted  in  this  region.  Unsteady  features  in  the 
flow  field  are  so  prevalent  that  it  is  unlikely  that  flow  stationarity  is  achieved  during  the  EFP’s  flight  time.  This  issue 
further  complicates  issues  associated  with  aero  stability.  We  have  also  investigated,  in  some  detail,  the  temperature 
distribution  on  the  exterior  and  interior  surfaces  of  the  EFP.  Purely  from  the  standpoint  of  aerodynamics,  areas  of 
severe  temperature  gradient  and  possible  thermal  shock  have  been  revealed.  We  have  also  identified  temperature 
regions  in  the  flow  field  that  may  support  the  existence  of  ionization  reactions  that  may  further  impact  the  dynamics 
of  this  system.  Given  the  turbulent  nature  of  the  flow  field,  ionization,  a  form  of  chemical  reaction,  should  be 
investigated.  We  plan  to  investigate  these  effects  in  future  work  associated  with  this  project.  The  present  studies  are 
limited  to  fixed  geometric  body  configurations;  however,  we  plan  on  implementing  new  algorithms  to  incorporate 
time-dependent  deformation  of  the  EFP  body  in  the  near  future. 
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